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A significant problem facing next-generation quantum technologies is how to generate and ma¬ 
nipulate macroscopic entanglement in light and matter systems. Here we report a new regime of 
dynamical light-matter behavior in which a giant, system-wide entanglement is generated by varying 
the light-matter coupling at intermediate velocities. This enhancement is far larger, broader-ranged, 
and more experimentally accessible, than that occurring near the quantum phase transition of the 
same model under adiabatic conditions. By appropriate choices of the coupling within this interme¬ 
diate regime, the enhanced entanglement can be made to spread system-wide or to reside in each 
subsystem separately. 
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I. INTRODUCTION 

Many-body quantum dynamics lie at the core of many 
natural phenomena and proposed quantum technologies, 
including information processing through schemes such 
as adiabatic quantum computing [T]. Achieving the con¬ 
trollable generation and manipulation of entanglement 
over many qubits is a key challenge, while doing so in 
light-matter systems is highly desirable for optoelectronic 
implementations. The ground state, and hence entan¬ 
glement, of a quantum system can be varied in a con¬ 
trolled way through adiabatic perturbations, though this 
is in principle an infinitely slow process. Quantum Phase 
Transitions (QPTs) can provide a naturally occurring en¬ 
tangled state and it is known that the entanglement can 
be enhanced at the critical point [2] . Recent studies have 
focused on time-dependent perturbations around QPTs 
that are either very slow (adiabatic) [3] or very fast (sud¬ 
den quench) [4]; or small dynamic oscillations around a 
phase space region [5]; or static coupling after a sudden 
quench im. 

Here we consider, by contrast, the regime of interme¬ 
diate perturbation velocities that has so far been over¬ 
looked. We consider the experimentally realized light- 
matter system of the Dicke Model (DM) [10], which 
has been realized in a variety of systems (e.g. circuit 
QED m and cold atom settings m USD- We un¬ 
cover a level of quantum complexity that is far richer 
than either the adiabatic or fast-quench regimes. The 
system-wide entanglement is dramatically enhanced over 
the static or adiabatic QPT values. Our results extend 
current understanding of coupled light-matter systems 
beyond the equilibrium ground state EHUHn], and also 
they also stand apart from more recent studies of out-of¬ 
equilibrium critical behavior EiiHailli. Moreover, our 
fully quantum analysis covers all dynamical regimes from 
very slow adiabatic through to sudden quench, capturing 
at each stage the emergent non-linear self interactions 
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and correlations within and between each subsystem. 
Our calculations employ the DM Hamiltonian m-- 


H = er 


■ 2Mlj (flt 

Viv 


( 1 ) 


where N is the number of matter qubits, the opera¬ 
tors Ji = ^ denote collective operators of the 

qubits, and operator (d) is the creation (annihilation) 
operator of the radiation field. In the thermodynamic 
limit {N oc), the critical value of the light-matter cou¬ 
pling parameter A is Ac = /2 while its finite-A" equiv¬ 

alent is slightly different [mus]. We treat the Hamilto¬ 
nian exactly using a large basis set |20] and avoid com¬ 
mon simplifications such as rotating-wave or semiclassi- 
cal approximations. The total system evolves unitarily 
under Hamiltonian H, starting at t = 0 from the instan¬ 
taneous ground state at A = 0. The light-matter coupling 
parameter is characterized by an annealing velocity (AV) 
V under a linear ramping A(t) = vt. More complicated 
time-dependencies can be treated but further complicate 
the understanding. The interval of interest in this paper 
is A G [0,2], meaning that the driving passes across a 
broad range of coupling strengths below and above the 
QPT. While the evolution of the total system S is de¬ 
scribed by a pure state |4^(t)), any subsystem A is de¬ 
scribed by a density matrix pA defined as the trace with 
respect to the other degrees of freedom not present in 
A: pA{t) = trs'-A (l^(^)) (^(^)l)- Because of the global 
unitary condition and Schmidt decomposition, both ra¬ 
diation and matter have the same value of entropy Sn, 
and it provides a quantitative measure of the degree of 
entanglement between them |2T] . 

This work is structured as follows: In section El we 
present the main result of our work, namely an enhanced 
dynamical light-matter regime at intermediate annealing 
velocities. It is followed, in section [in| by a deeper theo¬ 
retical explanation of our findings in terms of a dynamical 
symmetry-breaking and effective non-linear interactions. 
Then, in section El we establish the prevalence of our 
results for a wide range of system sizes and even if the 
system is submitted to a dissipating environment. Fi- 
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nally, we provide some concluding remarks in section [V| 


II. ENHANCED DYNAMICAL LIGHT-MATTER 
ENTANGLEMENT 


Figure summarizes our main finding: A novel, 
dynamical light-matter regime with greatly enhanced 
system-wide properties including entanglement, when 
the light-matter coupling is driven at intermediate veloc¬ 
ities. In results of Fig. [^and hereinafter, we are assum¬ 
ing resonant energies between qubits and held mode, i.e., 
u; = e in Eq. The peak entanglement value (purple) is 
far larger than the known equilibrium critical maximum 
[2], i.e. much larger than what can be achieved under 
adiabatic conditions. Also, this enhanced entanglement 
extends over a far broader range, well across the A > Ac 
region. As the annealing velocity increases, the critical 
onset point of light-matter entanglement is pushed to¬ 
ward larger A values and is no longer represented by a 
sharp peak, but instead a wavy plateau. At much higher 
velocities beyond the giant entanglement regime, A varies 
so fast that a sudden quench condition is achieved. Now 
the system is not quick enough to respond to the light- 
matter coupling, at least not in the A G [0, 2] interval of 
Fig. 0 

Besides exhibiting greatly enhanced values of light- 
matter entanglement, the intermediate regime has also 
some crucial practical advantages. Dotted lines in Fig. 

show how different are the time scales for optimal 
maximum light-matter entanglement in the intermediate 
regime as compared to the adiabatic one. In the adiabatic 
regime, the optimal value of light-matter entanglement is 
at the critical point Ac- By contrast, in the intermediate 
regime the optimal values are achieved well inside the 
A > Ac phase and require evolution times that are just 
a tiny fraction of the time needed to reach the critical 
point in the adiabatic regime. In any realistic implemen¬ 
tation, the system losses quantum information towards 
the environment. If these open system effects are not 
negligible, the intermediate regime is the only viable way 
to achieve optimal light-matter entanglement, before the 
dissipation of the environment becomes relevant. 


HI. DYNAMICAL SYMMETRY-BREAKING 
AND EFFECTIVE NON-LINEAR 
INTERACTIONS 


We now develop a deeper theoretical understanding of 
the results in Fig. by analyzing the underlying quan¬ 
tum state in the three main dynamical regimes, as il¬ 
lustrated in Fig. We start by rewriting the Dicke 
Hamiltonian exactly as 


4A2 

H = cob^b - + 
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( 2 ) 



Annealing parameter, A(t) 


FIG. 1. (color online) Dynamic evolution of the Von- 
Neumann entropy Sn (time varies from left to right) for the 
Dicke Model (DM) with iV = 81 qubits. For simplicity we set 
e = Ct; = 1 in Eq. [^for this hgure, as well as for Figs, [^and 
1^ The velocity range spans all regimes: from adiabatic (bot¬ 
tom) to sudden quench (top). Roman numerals mark where 
the held (Wigner) and matter (Agarwal-Wigner) distributions 
are depicted in Fig. Dotted lines represent instants of equal 
equal time t passed after initial condition (from up to down: 
t = 2,4,8,16). In between those curves there is a region of 
novel dynamical light-matter behavior with greatly enhanced 
entanglement (purple) as compared to the QPT. It arises for 
intermediate velocities and lies well inside the coupling regime 
for the conventional ordered phase (A > Ac = 0.5). The QPT 
corresponds toA^O.Sasr'^O and hence as log v ^ —oo 
(i.e. it tends toward A = 0.5 on the horizontal axis of the 
diagram). All dimensional quantities throughout this work 
are expressed in units so that energies in Eq. ^are cc = e = 1. 
Inset: a schematic representation of the DM which mimics 
various experimental realizations. 


where b = a In the A > Ac range, the last term 

becomes less and less relevant and the Dicke Hamilto¬ 
nian can be seen as a radiation mode that feels a dis¬ 
placed harmonic potential whose values depend on the 
eigenstate \mx) of in which the matter system sits [4]. 
Specifically, if A ^ Ac, then 



p + [x — 


2A 

— 1 =^ 
oVN 



\mx) {mx \ , 


(3) 

where we have used the quadrature operators of the radi¬ 
ation mode. The different confining potentials depending 
on the eigenvalue of are depicted by different parabo¬ 
lae in Eig. Importantly, the energy potential is sym¬ 
metrical with respect to a change in sign in which is 
a source of degeneracy. In addition, as \mx\ gets bigger, 
the minimum value of the harmonic potential becomes 
lower. The ground state of this approximate Hamilto¬ 
nian is any superposition of the form 


IV’o) = cos6>|A/2)„^^ |-/?)+e*‘^sin6>|-A/2)„^^ |/3), (4) 
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i.e. it corresponds to the two minimum parabolae. Hence 
the symmetry of the ground state is spontaneously bro¬ 
ken. The field state |/3) is a coherent state with P = 
As parity is preserved during the ramping, and 

ujy N 

adiabatic evolution keeps the energy in the lowest possi¬ 
ble value, the projection of |'0o) onto the even parity sub¬ 
space {0 = 7r/4 and (^ = 0 in Eq. is the state achieved 
in this regime. Both the symmetry breaking and the adi¬ 
abatic asymptotic value of entropy {Sn = log 2) in each 
subsystem can be explained by this double-well, since two 
coherent states are needed to describe each system. 


state in the intermediate regime. Other forms of dynam¬ 
ically enhanced non-adiabatic entanglement generation 
are possible, including the dynamical evolution preceded 
by a sudden quench [22] . 



FIG. 2. (color online) Schematic indicates energies of the 
three main outcomes when the light-matter coupling is in¬ 
creased at different velocities. Horizontal axis is the field posi¬ 
tion quadrature, both for sketches of the field Wigner function 
and for the harmonic confining potentials (parabolae). The 
potential felt by the field depends on the eigenvalue rrix of the 
matter operator Jx- If driving is very slow, the system stays 
cool in a symmetry broken ground state (blue curve) where 
each subsystem has entropy log 2 (low matter-light entangle¬ 
ment). If driving is very fast (sudden quench), the system 
still is in the initial A = 0 state (red curve) and light-matter 
entanglement is zero. For intermediate velocities, the system 
gets heated such that it occupies a complicated superposition 
of all the Jx eigenvalues, and hence gets disturbed by all these 
potentials (yellow curve). 

The multiple potential picture in Fig. and Eq. 
is still valid for the intermediate and sudden quench 
regimes. As the AV increases, the process generates a 
relative heating with respect to the ground state. Eor 
high enough AV (sudden quench), the system stays es¬ 
sentially in its starting condition, and the heating is just 
the consequence of the initial state being very different 
from the instantaneous ground state. Despite this sudden 
heating being very high, the simplicity of the initial state 
leads to no matter-light entanglement (i.e. one coherent 
state describes each subsystem). In the novel intermedi¬ 
ate regime, by contrast, all confining potentials simulta¬ 
neously perturb the system. A complicated superposition 
of non-trivial states for each parabola is generated, with 
complex and chaotic features. The distribution of prob¬ 
abilities shows a complex distribution across Jx eigen¬ 
values, leading to an entropy of each subsystem that is 
significantly higher in the intermediate regime than in the 
other two regimes. Due to its complex nature, there is 
no simple way to describe the structure of the dynamical 


FIG. 3. (color online) (a,b) Field Wigner distributions at 
the same instants as the Roman numerals marked in Fig. 

(I) becomes (a) and (H) becomes (b). Horizontal scales are 
for p momentum quadrature, while vertical ones are for the 
X position quadrature. Both figures have same x — p scale 
but the color scale is very different, —x part of the distri¬ 
bution not shown because it is symmetrical due to parity: 
W{x,p) = W{—x,—p). (c,d) Respective Agarwal-Wigner 

distributions of the qubits in the Bloch sphere. Negative re¬ 
gions and fragmentation (yellow curve of Fig. are a char¬ 
acteristic of the giant entanglement regime. Opposite hemi¬ 
spheres are not shown, but can be inferred from the symmetry 
relation W{0, (j) tt) = W{0, 0) due to conserved parity. 

The distinct behavior in the intermediate regime can 
also be interpreted in terms of the realization of non¬ 
linear self-interactions in each subsystem Haim, with 
qubits subjected to One-Axis Spin Squeezing with a 
Transverse Field [23l|24]. As the light-matter coupling 
increases in time from zero into the A > Ac range, the 
system moves out of its frozen initial state. The interac¬ 
tion term in the Dicke Hamiltonian begins to dominate 
and each subsystem works as a mediator of the other’s 
self-interaction. Matter-light interaction is in principle 
linear and generates little entanglement in the adiabatic 
limit, since the system has enough time to continuously 
stabilize in response to the perturbation. The interme¬ 
diate regime is the only one in which non-linearities can 
develop significantly. High values of squeezing in both 
subsystems are then generated [20]. However, this pro¬ 
cess does not last indefinitely: each subsystem starts to 
retain quantum information of the other, which is the mo¬ 
ment when the entropy grows. The effective interaction 
is broken, leaving the qubits and radiation mode entan¬ 
gled with each other but not within themselves. Despite 
being a single radiation mode, the field acts as a reservoir 
that dissipates the quantum correlations present in the 
squeezed states of each subsystem. The quantum corre¬ 
lations developing in the field system can be represented 
by the Wigner quasi-probability distribution in Figs. ia) 
and (b). See appendix [A| for definitions of these distribu- 
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tions. This shows that the distribution becomes highly 
fragmented yet retains some order, reflecting the com¬ 
plexity of the light subsystem in the intermediate regime. 
Round-tailed interference patterns as in Fig. [^have been 
obtained in light with a non-linear Kerr-like interaction 
following a Fokker-Planck equation [25]. This confirms 
that the field experiences an effective non-linear interac¬ 
tion. Similar signatures of complexity arise in the mat¬ 
ter subsystem, specifically the matter density matrices, 
as shown in the spherical Agarwal-Wigner functions in 
Figs, [^c) and (d) [26] . 



FIG. 4. (color online) Dynamical phase diagram showing the 
new enhanced entanglement regime (shaded red) defined by 
instances where Sn > log 2 since Sn ^ log 2 is the asymptotic 
value for the equilibrium ordered phase. For visual clarity in 
this figure, we take Sn > log 2 -h 0.05. Below a minimum 
annealing velocity Umin the system’s behavior is adiabatic, 
while for v > Umax it corresponds to sudden quench. The 
adiabatic boundary depends on system size as Umin oc N~^. 
The sudden quench boundary is size-independent, but de¬ 
pends on the value of Xd that is reached during the driving: 
Umax oc {Xd — Xc)^^^. See appendix [ b| for more details. 


IV. PERSISTENCE OF RESULTS AGAINST 
DISSIPATION AND DIFFERENT SYSTEM SIZES 

Figure [^demonstrates how the range of velocities that 
classify as ‘intermediate’ actually increases with increas¬ 
ing number of qubits V, meaning that the enhanced en¬ 
tanglement regime (EER) begins to dominate the space 
of behaviors as opposed to becoming a small niche. The 
EER can be imagined as lying between a lower bound 
AV Umin which marks the adiabatic evolution, and an 
upper bound one Umax defining the AV at which the sud¬ 
den quench approximation starts to be valid. Specifically, 
Fig. 11 shows the dynamical phase diagram of the inter¬ 
mediate regime in which the giant entanglement occurs, 
including its scaling behavior (dependence on system size 
N). The adiabatic evolution is more difficult to achieve 


as the number of atoms increases. The other main vari¬ 
able is the value of Xd reached by the annealing. The 
sudden quench condition requires higher AVs as this Xd 
gets bigger. The oscillatory behavior near the adiabatic 
regime has been smoothed out in order to make the phase 
boundary visually clearer. More details on how this dia¬ 
gram was obtained are given in appendix [B| 

In section [11] we used Von Neumann entropy S'at as 
a measure of light-matter entanglement in the case that 
the entire system is pure. However, in order to have an 
idea of the decoherence effects of a leaky cavity we have 
to analyze the DM as an open system, and then, Sn 
no longer a good entanglement witness. Instead, we use 
quantum negativity, whose non-zero value is a sufficient 
condition for bipartite entanglement in the open light- 
matter system m- Quantum negativity is defined as, 

115^"^ II — 1 

-, (5) 


where is the partial transpose of p with respect to 
the matter subsystem, and A = tr 


{VTi} i 


is the 


trace norm. However, one may wonder if switching from 



Annealing parameter, A(t) 


FIG. 5. (color online) Dynamical profile analogous to Fig. [^ 
but now for logarithmic negativity log 2 (2A/’+ 1). There is 
strong resemblance between both profiles. This justifies the 
use oi Sn as the entanglement witness for the closed system 
case while changing to log 2 {2J\f -hi) when the system is open. 
The inset shows the dynamical evolution of logarithmic neg¬ 
ativity for a near-adiabatic ramping velocity log 2 (u) = —8.95 
(dashed lines) and an intermediate one log 2 ('u) = —1.55 (solid 
lines). Different line colors represent different system sizes. 
This inset shows that even with relative small system sizes, 
log 2 (2A/’ -h 1) has qualitatively similar behaviors for the same 
u, so that conclusions about the robustness of the dynamic 
light-matter entanglement against decoherence can be extrap¬ 
olated to bigger N. 


one form of entanglement measure to the other has any 
justification. Figure [^ shows a dynamical profile of the 
ramping process analogous to that of Fig. [^ (same pa¬ 
rameters, unitary evolution, and system size N = 81), 
but with logarithmic negativity log 2 (2A/'+ 1). The re¬ 
semblance with the Sn graphic is quite apparent, and 
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the EER is again clearly noticeable. This should be no 
surprise as both measures virtually codified the same in¬ 
formation. If total p is a pure state, and {pi} is the 
spectrum of pq in that case (which is the same as the 
spectrum of because of Schmidt decomposition), both 
negativity JV and Sn can be written in terms of that 
spectrum. The first one would be 2Af + 1 = Xli j -\/PiPj ^ 
while the second one is Sn = —"^iPi^ogpi. Therefore, 
Fig. [5] confirms that both entanglement measures are 
well connected and suitable for witnessing light-matter 
entanglement. The reason why Sn was preferred in sec¬ 
tion [n] is because its wider usage and well established 
connection to other quantum information concepts [21]. 



FIG. 6. (color online) Effect of cavity losses on light- 
matter entanglement generation, witnessed by quantum log¬ 
arithmic negativity. Cavity’s field decay rate is k (different 
line style). Black curves depict giant entanglement regime 
log 2 ('r’) = —1.58, while red curves represent adiabatic one 
log 2 (n) = —8.96. Main figure is for iV = 5 system’s size 
while inset represents analogous N = 11 results. Dynamical 
entanglement regime is bigger and more robust against open 
system losses as N increases. In opposition, near-adiabatic 
entanglement vanishes even for weak losses since annealing 
time is too long as compared to decoherence time. Zero tem¬ 
perature has been assumed and k. is measured in units of field 
free frequency uj. Finite temperature cases have very similar 
tendencies. 

It should be stressed that the giant entanglement 
regime, central to present work, is accessible under cur¬ 
rent experiments very similar to the sketch in inset of Fig. 

In particular, Klinder et al. have been able 
to ramp a cold atom Bose-Einstein condensate whose dy¬ 
namics are governed by the DM, and with annealing ve¬ 
locities corresponding very well to the giant entanglement 
regime predicted by our analysis m- Moreover, cavity 
decay rates, which measure the importance of losses, are 
an order of magnitude lower than the main evolution en¬ 
ergy cj, and can be well simulated by open cavity mas¬ 
ter equations [13]. Our analysis including losses predicts 
the survival of giant light-matter entanglement in those 
experiments - indeed it occurs even when open system 
losses are present, as shown in Eig. We use quantum 
negativity as the light-matter entanglement measure for 
the open system case m since Sn will no longer be a 


good bipartite entanglement measure. Details of the def¬ 
inition of quantum negativity and the open system solu¬ 
tion, are given in the appendix Since the light-matter 
entanglement increases as N increases, we predict that 
experimental system sizes N ^ 10^ will develop signif¬ 
icantly more robust light-matter entanglement against 
losses than for small systems. By contrast, entanglement 
generation in the well-studied near-adiabatic regime can¬ 
not be achieved under current experimental setups, since 
decoherence times are much shorter than the annealing 
times required by this regime. 

We have then established the persistence of the EER 
for a wide range of system sizes and even under the ef¬ 
fect of dissipation effects. A remnant concern could be 
that the system sizes N accessible to open system numer¬ 
ical solution in the present work are well smaller than the 
ones examined in the pure case. This could cast doubt on 
whether the conclusions brought from this small size re¬ 
sult have any general validity. The inset of figure [^ shows 
that despite even at N = 11, there are good signatures 
of the dynamical phase with giant light-matter entangle¬ 
ment, and that it only gets stronger as N increases. This 
means that, provided the open system’s parameters like 
K and n remain the same, any deleterious effect caused 
by the leakage of quantum information to the environ¬ 
ment is only weaker with sizes of the order A" = 81 than 
the ones numerical accessible. 


V. CONCLUSIONS 

We have unraveled a previously unnoticed dynamical 
regime where enhanced light-matter entanglement arises 
at intermediate ramping velocities. We have provided 
theoretical insights of this phenomenon by means of a 
dynamical symmetry breaking and effective non-linear 
self-interactions. We have further argued that our main 
results are still present under more realistic settings, as 
those already accessible by current experiments, with big¬ 
ger system sizes and a dissipating environment. Einally, 
we note that by moving around the parameter space in 
time in Eig. the enhanced squeezing and entanglement 
can be altered within the matter and light subsystems 
separately, and then transferred by means of the light- 
matter coupling. Potential applications include high pre¬ 
cision quantum metrology and a range of quantum infor¬ 
mation processing technologies [28H3Q] . 
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Appendix A: Subsystem quasi-probabilities 
distributions 

In order to provide a visual depiction of both matter 
and field subsystems in the main text, we employed phase 
space representations of the density matrix by means of 
Wigner functions. The matter states were represented by 
the Agarwal-Wigner function (AWF), which is a Bloch 
sphere representation of the subsystem’s density matrix 

Pq [2Sj 5 


N I 

Wg{e,<l)) = (Al) 

Z=0 m=—l 

where Yi^rn are the spherical harmonics; and terms Ti^rn = 
trjpqTi.m} are expected values of the multipole opera¬ 
tor, 

fi,m= Yl 

M,M' = -j 

(-M i 

where j = A^/2, and ( ^ ^/ ) is the Wigner 3j 

symbol. 

The respective Wigner function for the field density 
matrix is HEI], 

oo 

Wb (a, Pb) = Y (-1)" («) PbD {a) \n) , (A3) 

n=0 


where D {a) = e^" " is the displacement operator, 

and Of G C. The displacement parameter a can be 
expressed in terms of the field’s position {x) and mo¬ 
mentum (p) quadratures as its real and imaginary part, 

\f20L = X + Ip. 


Appendix B: Numerical evidence of phase diagram 

This section provides a detailed explanation on how we 
obtained the phase diagram of the EER depicted in Eig. 
[^of main text. The EER was defined as the dynamical 
region in Eig. □of main text where entropy is bigger than 
log 2 + 0.05 well inside the ordered A > Ac phase. This re¬ 
gion can be seen as bounded by a maximum AV 'L’max over 
which the sudden quench approximation is valid, and a 
minimum AV 'L’min below which the adiabatic condition 
is fulfilled. 

Erom Eig. of main text, it is clear that 'L’min does 
not depend on the value Xd of the annealing parameter 
reached (the lower boundary forms an horizontal line in 
the figure). This is expected as the ground state in all 
the ordered phase has an asymptotic value of S'at -> log 2, 
and the adiabatic condition should only depend on the 



Log. of AV: \og^(v)or \og^(Nv) 


FIG. 7. (color online) Numerical evidence for constructing 
the dynamical phase diagram of Fig. of main text. Solid 
(dashed) curves depict the value of aSat as a function of log 2 ('c) 
(log 2 {Nv)) when A = 2.0 for different system sizes. Point 
'Cmax (iV'Cmin) is obtained in a size-independent way where 
curves touch the red horizontal line, which is the defining 
value for the dynamical phase boundary. The inset shows the 
dependence of point 'Umax on the value of Xd reached, and its 
fitting to the Umax oc (A —Ac)^^^ dependence (red dashed line). 
In opposition, the point Nvmin is independent of the value of 
Xd reached. In principle, solid and dashed curves represent the 
same data, but dashed curves have been smoothed in order 
to establish a cleaner boundary for Vmin- 


system size N. The scaling 'L’min oc N~^ comes from a 
well established relation of the minimal energy gap A 
at the critical threshold 0 US], and it is confirmed in 
the dashed curves of Fig. where using this kind of 
scaling of the AV makes the lower boundary of the EER 
size independent.The upper bound, 'C’max depends on Xd 
and the inset of Eig. shows that the fitting 'Umax oc 
{Xd — Ac)^/^ is a good approximation. On the other, it 
is clear that 'Umax does not depend on system size, as all 
three solid curves of Eig. [preach the red line at the same 
point. The relations Umax oc {Xd — Xc)^^^ and Umin oc N~^ 
with the addition of one particular point (as the ones 
obtained in Eig. is all that is needed to produce Eig. 
3 in main text. 


Appendix C: Open System Evolution 


The composed DM light-matter density matrix p{t) 
evolves under a unitary part generated by H (see Eq. 1 
of main text) and a dissipative part due to a field lossy 
cavity. The evolution of p{t) is modeled by a Master 
Equation in Lindblad form [32] . 


—p = —L 

di 


H,p 


2k {n 1)7) (p; a) + 2KnV (p; d^) , 

(Cl) 
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where k is the damping rate of the cavity, and n is the 
thermal mean photon number. Also, for any operator A, 
the Lindblad superoperator V is defined as ({,} denotes 
anti-commutator), 

v[p-,A) = ApA^-^{A^A,p]. (C2) 

The general initial condition (when A(0) = 0) consist of 
an unentangled light-matter state, 

p-/3a'*'a 

^(O) = |-W/2).(-»/2|.0 ^^.j^.,.,.^ , (C3) 

where e~^ = n/ (n-1-1), and /3 = 1/T is the inverse 
temperature in natural units. This condition means zero 


excitations in the matter state and field thermalization. 
The field state is in fact the steady equilibrium state for 
the respective terms in Eq. |C1[ 

We have obtained exact numerical solutions of Eq. |C1[ 
either in the pure DM Hilbert state when /i: = 0, or in 
the space of density matrices otherwise, using the DM 
Hamiltonian or the Liouvillian superoperator as the re¬ 
spective dynamical generator [32]. Since the dimension 
of the numerical evolution vector space is squared as soon 
as the unitary condition is broken, only relatively small 
values of N can be investigated under open conditions. 
As the field Hilbert space is infinite dimensional, we have 
truncated its Eock basis {|n)} up to a number where nu¬ 
merical results converge. 
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